!
!
! Roel Matthysen
! TWS: Huistaak 6
!
!
! Gebaseerd op de code gebruikt tijdens de oefenzitting
!

program multest
	! module zoals opgesteld in de oefenzitting
    use matrixop
    implicit none
    
    !--------------------------------------------------------------------------
    ! Abstract interfaces
    !
    ! NOTE: this simplifies the timings, but is out of the scope of the course.
    !--------------------------------------------------------------------------

    abstract interface
        subroutine a_maal_b_interface(a, b, c)
            import dp
            real(kind=dp), dimension(:,:), intent(in)  :: a, b
            real(kind=dp), dimension(:,:), intent(out) :: c
        end subroutine a_maal_b_interface
        subroutine a_maal_b_blocks_interface(a,b,c,blockize)
            import dp
            real(kind=dp), dimension(:,:), intent(in)  :: a, b
            real(kind=dp), dimension(:,:), intent(out) :: c
            integer, intent(in) :: blockize
        end subroutine a_maal_b_blocks_interface
    end interface
    
    !--------------------------------------------------------------------------
    ! Main timing program
    !--------------------------------------------------------------------------

    integer :: k, N
    real :: flops
    integer, dimension(:), allocatable :: seed
    real(kind=dp), dimension(:,:), allocatable :: a, b, c
    real(kind=dp), dimension(:,:), allocatable :: c_matmul
    

    ! Make sure we use the same pseudo-random numbers each time by initializing
    ! the seed to a certain value.
    call random_seed(size=k)
    allocate(seed(k))
    seed = N
    call random_seed(put=seed)


    ! Start the timings
    print *, ""
    write(unit=*, fmt="(A)") "TIMING RESULTS:"
    
    ! 1. Three nested loops
    call do_timing( "JKI", a_maal_b_jki )
    
    ! 2. Two nested loops with dot_product and explicit transpose of matrix A
    call do_timing( "JI TP DOT_PRODUCT", a_maal_b_transp_ji_dot_product )
    
    ! 3. Using BLAS
    call do_timing( "BLAS DGEMM", a_maal_b_blas )
    
    ! 4. In blocks
    call do_timing( "IN BLOCKS", method_blocks=a_maal_b_blocks )
    
    ! 5. Intrinsic matmul function
    call do_timing( "MATMUL", a_maal_b_matmul )
    
contains

    subroutine do_timing( name, method, method_blocks )
    
        character(len=*), intent(in) :: name
        procedure(a_maal_b_interface), optional :: method
        procedure(a_maal_b_blocks_interface), optional :: method_blocks
        
        character(len=len_trim(name)+4) :: fnameout
        integer, parameter :: nb_repeat = 3
        integer :: i_repeat, i_N
        real(kind=dp) :: mynorm
        real :: t1, t2, t(3)
        ! Opstellen van de testwaarden voor N
        integer :: Ns(8) = (/ (2**(i_N+2), i_N=1,8)  /)
        
        write(0,'(a)') name
        
        
        
        WRITE(fnameout,'(a)') name//".out"
		OPEN(33,file=fnameout)	
        
        do i_N = 1,size(Ns,1)
        
            ! Matrix multiplication is an O(N^3) operation if carried out naively, because
    		! for each of the N^2 elements of the result, we perform N multiplications and
    		! N-1 additions.  Together, this is N^2*(N+N-1) flops.
        	N = Ns(i_N)
            flops = real(N)*N*(2*N-1)
    		write(unit=*, fmt="(A, I0, A)") &
        	"  Computational cost: ", int(flops/1.0e6), " Mflop."
			! Allocate the matrices itself and one reference matrix
			allocate(a(N,N), b(N,N), c(N,N))
			call random_number(a)
			call random_number(b)
			call a_maal_b_matmul(a,b,c_matmul) ! Reference value
			write(unit=33, fmt="(i6)", advance="no") & 
				N
			! Do the timing (3 maal)
			do i_repeat = 1,nb_repeat
				if( present(method) ) then
					call cpu_time(t1)
						call method( a, b, c )
					call cpu_time(t2)
				else
					call cpu_time(t1)
					! Optimale blocksize wordt gebruikt, anders blocksize = N
						call method_blocks( a, b, c, min(512,N))
					call cpu_time(t2)
				end if
				t(i_repeat) = (t2-t1)
				if( t(i_repeat) /= 0 ) then
					write(unit=33, fmt="(ES9.2)", advance="no") &
						(flops/t(i_repeat))/(1000*1000)
				else
					write(unit=33, fmt="(A9)", advance="no") &
						"Inf"
				end if
			enddo
	
			! Output the result
	
			write(unit=33, fmt="(A)") " "            
			
			! Clean up
			deallocate(a, b, c)
    	
		enddo
        
    end subroutine do_timing

end program multest
